%% Setting Directory & Paths
clear
cd /yourpath/FinalReplicationPackage/

%% Define Parameters
b_call = 0.0949;
b_weather = -0.0724;
gamma = 0.484;
%gamma = 0.692;
plot_figure2 = "yes"

N_c0 = 350;
N_c1 = 350;
N_c2 = 350;
N_c12 = 350;

sargan_p_array = [];
d_past_coef_array = [];
R_present_coef_array = [];
P_present_coef_array = [];

%% Call Indicator
P_c1_1 = ones(1,N_c1);
P_c2_1 = zeros(1,N_c2);
P_c12_1 = ones(1,N_c12);
P_c0_1 = zeros(1,N_c0);

    P_i1 = [P_c1_1, P_c2_1, P_c12_1, P_c0_1];

P_c1_2 = zeros(1,N_c1);
P_c2_2 = ones(1,N_c2);
P_c12_2 = ones(1,N_c12);
P_c0_2 = zeros(1,N_c0);

    P_i2 = [P_c1_2, P_c2_2, P_c12_2, P_c0_2];

%% Define alpha = 1, i.e. no motivation based persistence
alpha = 1

%% Define delta, i.e. cyclical change in baseline motivation
delta_vec = [1]

frac_B = []
frac_b = []
frac_delta = []

eta_vec = [0.05]

rand('seed', 1647853);

for eta = eta_vec
    for delta = delta_vec
    q = 0.5  
    binary = [delta, -delta];
    weights = [1-q, q]; 
    
        delta_i = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);
        eta_i = eta*((delta_i - 1)/2)
        eta_j = -eta*((delta_i + 1)/2)
    
    %% 10000 loops
    %for n = 1:10000
    
        %n;
    
    %% Adverse Weather Indicator
    binary = [0, 1];
    weights = [1-.204, .204]; % (20.4) 23.3% of Adverse Weather
    
        R_i1 = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);
        R_i2 = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);
        R_i3 = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);
        R_i4 = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);
    
    % Draw epsilons
    e_i1 = rand(1,N_c1+N_c2+N_c12+N_c0);
    
    % Selecting b_0, such that mean donation of c0 = 0.1275
    b_0_choice_array = 0:0.001:1;
    
        for i = 1:length(b_0_choice_array)
            b_0_choice(i,1) = b_0_choice_array(i);
            u_i1_c0 = b_0_choice_array(i)+(b_0_choice_array(i)*delta_i(1051:1400))+eta_i(1051:1400)+(b_weather*R_i1(1051:1400));
            mean_donation_control(i,1) = mean(u_i1_c0>e_i1(1051:1400));
        end
    
    diff = abs(mean_donation_control-0.1275);
    min_diff=min(abs(mean_donation_control-0.1275));
    b_0 = mean(b_0_choice(diff==min_diff));
    B = b_0;
    
    %% Simulating Period 1
    % Compute Utility
    u_i1 = B + B*delta_i + eta_i + b_call*P_i1 + b_weather*R_i1;
    
    % Extract d_i1
    d_i1 = u_i1>e_i1;
    
    % Compute Mean Donation rate
    mean_c1_p1 = mean(d_i1(1:350));
    mean_c2_p1 = mean(d_i1(351:700));
    mean_c12_p1 = mean(d_i1(701:1050));
    mean_c0_p1 = mean(d_i1(1051:1400));
    
    se_c1_p1 = std(d_i1(1:350))/sqrt(N_c1);
    se_c2_p1 = std(d_i1(351:700))/sqrt(N_c2);
    se_c12_p1 = std(d_i1(701:1050))/sqrt(N_c12);
    se_c0_p1 = std(d_i1(1051:1400))/sqrt(N_c0);
    
    
    %% Simulating Period 2
    % Compute Utility
    u_i2 = B + B*delta_i*-1 + eta_j + b_call*P_i2 + b_weather*R_i2 + alpha*gamma*d_i1;
    % Draw Epsilons
    e_i2 = rand(1,N_c1+N_c2+N_c12+N_c0);
    
    % Extract Donation Indicator
    d_i2 = u_i2>e_i2;
    
    % Donation Rate
    mean_c1_p2 = mean(d_i2(1:350));
    mean_c2_p2 = mean(d_i2(351:700));
    mean_c12_p2 = mean(d_i2(701:1050));
    mean_c0_p2 = mean(d_i2(1051:1400));
    
    % Standard Error
    se_c1_p2 = std(d_i2(1:350))/sqrt(N_c1);
    se_c2_p2 = std(d_i2(351:700))/sqrt(N_c2);
    se_c12_p2 = std(d_i2(701:1050))/sqrt(N_c12);
    se_c0_p2 = std(d_i2(1051:1400))/sqrt(N_c0);
    
    %% Simulating Period 3
    % Compute Utility
    u_i3 = B + B*delta_i + eta_i + b_weather*R_i3 + alpha*gamma*d_i2;
    
    % Draw epsilons
    e_i3 = rand(1,N_c1+N_c2+N_c12+N_c0);
    
    % Extract d_i3
    d_i3 = u_i3>e_i3;
    
    % Compute Mean Donation rate
    mean_c1_p3 = mean(d_i3(1:350));
    mean_c2_p3 = mean(d_i3(351:700));
    mean_c12_p3 = mean(d_i3(701:1050));
    mean_c0_p3 = mean(d_i3(1051:1400));
    
    % Compute Standard Error
    se_c1_p3 = std(d_i3(1:350))/sqrt(N_c1);
    se_c2_p3 = std(d_i3(351:700))/sqrt(N_c2);
    se_c12_p3 = std(d_i3(701:1050))/sqrt(N_c12);
    se_c0_p3 = std(d_i3(1051:1400))/sqrt(N_c0);
    
    %% Simulating Period 4
    % Compute Utility
    u_i4 = B + B*delta_i*-1 + eta_j + b_weather*R_i4 + alpha*gamma*d_i3;
    
    % Draw epsilons
    e_i4 = rand(1,N_c1+N_c2+N_c12+N_c0);
    
    % Extract d_i4
    d_i4 = u_i4>e_i4;
    
    % Compute Mean Donation rate
    mean_c1_p4 = mean(d_i4(1:350));
    mean_c2_p4 = mean(d_i4(351:700));
    mean_c12_p4 = mean(d_i4(701:1050));
    mean_c0_p4 = mean(d_i4(1051:1400));
    
    % Compute Standard Error
    se_c1_p4 = std(d_i4(1:350))/sqrt(N_c1);
    se_c2_p4 = std(d_i4(351:700))/sqrt(N_c2);
    se_c12_p4 = std(d_i4(701:1050))/sqrt(N_c12);
    se_c0_p4 = std(d_i4(1051:1400))/sqrt(N_c0);
    
    %% Tabulation of Data 
    df = table;
    
        df.d1=d_i1';
        df.d2=d_i2';
        df.d3=d_i3';
        df.d4=d_i4';
        df.e1=e_i1';
        df.e2=e_i2';
        df.e3=e_i3';
        df.e4=e_i4';
        df.u1=u_i1';
        df.u2=u_i2';
        df.u3=u_i3';
        df.u4=u_i4';
        df.P1=P_i1';
        df.P2=P_i2';
        df.R1=R_i1';
        df.R2=R_i2';
        df.R3=R_i3';
        df.R4=R_i4';
        
    %% Plot Labels 
    xpt_c0 = [1+0.1, 2+0.1, 3+0.1, 4+0.1];
    xpt_c2 = [1, 2, 3, 4];
    xpt_c12 = [1+0.05, 2+0.05, 3+0.05, 4+0.05];
    xpt_c1 = [1-0.05, 2-0.05, 3-0.05, 4-0.05];
    
    ypt_c1 = [mean_c1_p1, mean_c1_p2, mean_c1_p3, mean_c1_p4];
    ypt_c2 = [mean_c2_p1, mean_c2_p2, mean_c2_p3, mean_c2_p4];
    ypt_c12 = [mean_c12_p1, mean_c12_p2, mean_c12_p3, mean_c12_p4];
    ypt_c0 = [mean_c0_p1, mean_c0_p2, mean_c0_p3, mean_c0_p4];
    
    se_c1 = [se_c1_p1, se_c1_p2, se_c1_p3, se_c1_p4];
    se_c2 = [se_c2_p1, se_c2_p2, se_c2_p3, se_c2_p4];
    se_c12 = [se_c12_p1, se_c12_p2, se_c12_p3, se_c12_p4];
    se_c0 = [se_c0_p1, se_c0_p2, se_c0_p3, se_c0_p4];
    
    ypt_c1_ci_ub = ypt_c1+se_c1;
    ypt_c1_ci_lb = ypt_c1-se_c1;
    ypt_c2_ci_ub = ypt_c2+se_c2;
    ypt_c2_ci_lb = ypt_c2-se_c2;
    ypt_c12_ci_ub = ypt_c12+se_c12;
    ypt_c12_ci_lb = ypt_c12-se_c12;
    ypt_c0_ci_ub = ypt_c0+se_c0;
    ypt_c0_ci_lb = ypt_c0-se_c0;
    
    %% Plot Donation Rate conditional on donation in period 1 and 2
    % Compute Mean Donation rate conditional on donation in P1 and P2
    mean_yes_p1_p1 = mean(d_i1(d_i1 == 1));
    mean_no_p1_p1 = mean(d_i1(d_i1 == 0));
    mean_yes_p1_p2 = mean(d_i2(d_i1 == 1));
    mean_no_p1_p2 = mean(d_i2(d_i1 == 0));
    mean_yes_p1_p3 = mean(d_i3(d_i1 == 1));
    mean_no_p1_p3 = mean(d_i3(d_i1 == 0));
    mean_yes_p1_p4 = mean(d_i4(d_i1 == 1));
    mean_no_p1_p4 = mean(d_i4(d_i1 == 0));
    
    N_yes_p1 = sum(d_i1 == 1);
    N_no_p1 = sum(d_i1 == 0);
    
    mean_yes_p2_p1 = mean(d_i1(d_i2 == 1));
    mean_no_p2_p1 = mean(d_i1(d_i2 == 0));
    mean_yes_p2_p2 = mean(d_i2(d_i2 == 1));
    mean_no_p2_p2 = mean(d_i2(d_i2 == 0));
    mean_yes_p2_p3 = mean(d_i3(d_i2 == 1));
    mean_no_p2_p3 = mean(d_i3(d_i2 == 0));
    mean_yes_p2_p4 = mean(d_i4(d_i2 == 1));
    mean_no_p2_p4 = mean(d_i4(d_i2 == 0));
    
    N_yes_p2 = sum(d_i2 == 1);
    N_no_p2 = sum(d_i2 == 0);
    
    % Compute Standard Error
    se_yes_p1_p1 = std(d_i1(d_i1 == 1))/sqrt(N_yes_p1);
    se_no_p1_p1 = std(d_i1(d_i1 == 0))/sqrt(N_no_p1);
    se_yes_p1_p2 = std(d_i2(d_i1 == 1))/sqrt(N_yes_p1);
    se_no_p1_p2 = std(d_i2(d_i1 == 0))/sqrt(N_no_p1);
    se_yes_p1_p3 = std(d_i3(d_i1 == 1))/sqrt(N_yes_p1);
    se_no_p1_p3 = std(d_i3(d_i1 == 0))/sqrt(N_no_p1);
    se_yes_p1_p4 = std(d_i4(d_i1 == 1))/sqrt(N_yes_p1);
    se_no_p1_p4 = std(d_i4(d_i1 == 0))/sqrt(N_no_p1);
    
    se_yes_p2_p1 = std(d_i1(d_i2 == 1))/sqrt(N_yes_p2);
    se_no_p2_p1 = std(d_i1(d_i2 == 0))/sqrt(N_no_p2);
    se_yes_p2_p2 = std(d_i2(d_i2 == 1))/sqrt(N_yes_p2);
    se_no_p2_p2 = std(d_i2(d_i2 == 0))/sqrt(N_no_p2);
    se_yes_p2_p3 = std(d_i3(d_i2 == 1))/sqrt(N_yes_p2);
    se_no_p2_p3 = std(d_i3(d_i2 == 0))/sqrt(N_no_p2);
    se_yes_p2_p4 = std(d_i4(d_i2 == 1))/sqrt(N_yes_p2);
    se_no_p2_p4 = std(d_i4(d_i2 == 0))/sqrt(N_no_p2);
    
    
    xpt_yes_p1 = [1+0.1, 2+0.1, 3+0.1, 4+0.1];
    xpt_no_p1 = [1, 2, 3, 4];
    
    xpt_yes_p2 = [1+0.05, 2+0.05, 3+0.05, 4+0.05];
    xpt_no_p2 = [1-0.05, 2-0.05, 3-0.05, 4-0.05];
    
    ypt_yes_p1 = [mean_yes_p1_p1, mean_yes_p1_p2, mean_yes_p1_p3, mean_yes_p1_p4];
    ypt_no_p1 = [mean_no_p1_p1, mean_no_p1_p2, mean_no_p1_p3, mean_no_p1_p4];
    ypt_yes_p2 = [mean_yes_p2_p1, mean_yes_p2_p2, mean_yes_p2_p3, mean_yes_p2_p4];
    ypt_no_p2 = [mean_no_p2_p1, mean_no_p2_p2, mean_no_p2_p3, mean_no_p2_p4];
    
    se_yes_p1 = [se_yes_p1_p1, se_yes_p1_p2, se_yes_p1_p3, se_yes_p1_p4];
    se_no_p1 = [se_no_p1_p1, se_no_p1_p2, se_no_p1_p3, se_no_p1_p4];
    se_yes_p2 = [se_yes_p2_p1, se_yes_p2_p2, se_yes_p2_p3, se_yes_p2_p4];
    se_no_p2 = [se_no_p2_p1, se_no_p2_p2, se_no_p2_p3, se_no_p2_p4];
    
    ypt_yes_p1_ci_ub = ypt_yes_p1+se_yes_p1;
    ypt_yes_p1_ci_lb = ypt_yes_p1-se_yes_p1;
    ypt_no_p1_ci_ub = ypt_no_p1+se_no_p1;
    ypt_no_p1_ci_lb = ypt_no_p1-se_no_p1;
    ypt_yes_p2_ci_ub = ypt_yes_p2+se_yes_p2;
    ypt_yes_p2_ci_lb = ypt_yes_p2-se_yes_p2;
    ypt_no_p2_ci_ub = ypt_no_p2+se_no_p2;
    ypt_no_p2_ci_lb = ypt_no_p2-se_no_p2;
    
    if plot_figure2 == "yes"
    
        if delta == delta_vec(1) & eta == eta_vec(1)
    
            clf
            %f = tiledlayout(2,3)
            %nexttile
    
        else
    
            %nexttile
    
        end  
    
        set(gca,'fontsize', 12.5)
      
        hold on
        f = plot(xpt_yes_p1(2:4), ypt_yes_p1(2:4), '-o', 'LineWidth', 1.5, 'MarkerSize', 10, color='r');
        plot(xpt_no_p1(2:4), ypt_no_p1(2:4), '-o', 'LineWidth', 1.5, 'MarkerSize', 14, color='r', marker="X");
        plot(xpt_yes_p2(3:4), ypt_yes_p2(3:4), '--o', 'LineWidth', 1.5, 'MarkerSize', 10, color='b');
        plot(xpt_no_p2(3:4), ypt_no_p2(3:4), '--o', 'LineWidth', 1.5, 'MarkerSize', 14, color='b', marker="X");
    
        plot(xpt_yes_p2(1), ypt_yes_p2(1), ':o', 'LineWidth', 1.5, 'MarkerSize', 10, color='b');
        plot(xpt_no_p2(1), ypt_no_p2(1), '-o', 'LineWidth', 1.5, 'MarkerSize', 14, color='b', marker="X");
        plot(xpt_yes_p2(1), ypt_yes_p2(1), ':o', 'LineWidth', 1.5, 'MarkerSize', 10, color='b');
        plot(xpt_no_p2(1), ypt_no_p2(1), '-o', 'LineWidth', 1.5, 'MarkerSize', 14, color='b', marker="X");
           
        plot(xpt_yes_p1(2:4), ypt_yes_p1_ci_ub(2:4), "_", 'LineWidth', 1.5, 'MarkerSize', 10, color='r');
        plot(xpt_yes_p1(2:4), ypt_yes_p1_ci_lb(2:4), "_", 'LineWidth', 1.5, 'MarkerSize', 10, color='r');
        plot(xpt_no_p1(2:4), ypt_no_p1_ci_ub(2:4), "_", 'LineWidth', 1.5, 'MarkerSize', 10, color='r');
        plot(xpt_no_p1(2:4), ypt_no_p1_ci_lb(2:4), "_", 'LineWidth', 1.5, 'MarkerSize', 10, color='r');
        plot(xpt_yes_p2(1,3:4), ypt_yes_p2_ci_ub(1,3:4), "_", 'LineWidth', 1.5, 'MarkerSize', 10, color='b');
        plot(xpt_yes_p2(1,3:4), ypt_yes_p2_ci_lb(1,3:4), "_", 'LineWidth', 1.5, 'MarkerSize', 10, color='b');
        plot(xpt_no_p2(1,3:4), ypt_no_p2_ci_ub(1,3:4), "_", 'LineWidth', 1.5, 'MarkerSize', 10, color='b');
        plot(xpt_no_p2(1,3:4), ypt_no_p2_ci_lb(1,3:4), "_", 'LineWidth', 1.5, 'MarkerSize', 10, color='b');
    
        plot(xpt_yes_p2(1), ypt_yes_p2_ci_lb(1), "_", 'LineWidth', 1.5, 'MarkerSize', 10, color='b');
        plot(xpt_yes_p2(1), ypt_yes_p2_ci_ub(1), "_", 'LineWidth', 1.5, 'MarkerSize', 10, color='b');
        plot(xpt_no_p2(1), ypt_no_p2_ci_lb(1), "_", 'LineWidth', 1.5, 'MarkerSize', 10, color='b');
        plot(xpt_no_p2(1), ypt_no_p2_ci_ub(1), "_", 'LineWidth', 1.5, 'MarkerSize', 10, color='b');
        
        plot([xpt_no_p2(1),xpt_no_p2(1)],[ypt_no_p2_ci_lb(1),ypt_no_p2_ci_ub(1)], color='b');
        plot([xpt_no_p2(2),xpt_no_p2(2)],[ypt_no_p2_ci_lb(2),ypt_no_p2_ci_ub(2)], color='b');
        plot([xpt_no_p2(3),xpt_no_p2(3)],[ypt_no_p2_ci_lb(3),ypt_no_p2_ci_ub(3)], color='b');
        plot([xpt_no_p2(4),xpt_no_p2(4)],[ypt_no_p2_ci_lb(4),ypt_no_p2_ci_ub(4)], color='b');
        
        plot([xpt_yes_p1(1),xpt_yes_p1(1)],[ypt_yes_p1_ci_lb(1),ypt_yes_p1_ci_ub(1)], color='r');
        plot([xpt_yes_p1(2),xpt_yes_p1(2)],[ypt_yes_p1_ci_lb(2),ypt_yes_p1_ci_ub(2)], color='r');
        plot([xpt_yes_p1(3),xpt_yes_p1(3)],[ypt_yes_p1_ci_lb(3),ypt_yes_p1_ci_ub(3)], color='r');
        plot([xpt_yes_p1(4),xpt_yes_p1(4)],[ypt_yes_p1_ci_lb(4),ypt_yes_p1_ci_ub(4)], color='r');
        
        plot([xpt_yes_p2(1),xpt_yes_p2(1)],[ypt_yes_p2_ci_lb(1),ypt_yes_p2_ci_ub(1)], color='b');
        plot([xpt_yes_p2(2),xpt_yes_p2(2)],[ypt_yes_p2_ci_lb(2),ypt_yes_p2_ci_ub(2)], color='b');
        plot([xpt_yes_p2(3),xpt_yes_p2(3)],[ypt_yes_p2_ci_lb(3),ypt_yes_p2_ci_ub(3)], color='b');
        plot([xpt_yes_p2(4),xpt_yes_p2(4)],[ypt_yes_p2_ci_lb(4),ypt_yes_p2_ci_ub(4)], color='b');
        
        plot([xpt_no_p1(1),xpt_no_p1(1)],[ypt_no_p1_ci_lb(1),ypt_no_p1_ci_ub(1)], color='r');
        plot([xpt_no_p1(2),xpt_no_p1(2)],[ypt_no_p1_ci_lb(2),ypt_no_p1_ci_ub(2)], color='r');
        plot([xpt_no_p1(3),xpt_no_p1(3)],[ypt_no_p1_ci_lb(3),ypt_no_p1_ci_ub(3)], color='r');
        plot([xpt_no_p1(4),xpt_no_p1(4)],[ypt_no_p1_ci_lb(4),ypt_no_p1_ci_ub(4)], color='r');
        
        xlabel("Period");
        ylabel("Donation Rate");
        xticks([1,2,3,4]);
        xlim([0.9 4.1]);
        ylim([-0.1 1.2]);
    
        lgd_lab1 = sprintf('Donation in P1 (N=%.0f)', N_yes_p1)
        lgd_lab2 = sprintf('No donation in P1 (N=%.0f)', N_no_p1)
        lgd_lab3 = sprintf('Donation in P2 (N=%.0f)', N_yes_p2)
        lgd_lab4 = sprintf('No donation in P2 (N=%.0f)', N_no_p2)
    
        lgd = legend(lgd_lab1, lgd_lab2, lgd_lab3, lgd_lab4, 'Location', 'Northwest');
        title(['delta=',num2str(delta) ', eta=', num2str(eta) ', gamma=', num2str(gamma)]);
        
    
        %filename = sprintf('output/figures/conditional_alpha=%.1f-gamma=%.3f-q=%.2f.png', alpha, gamma, q);
        %saveas(f,filename);
        
        frac_delta = [frac_delta, delta]
        frac_B = [frac_B; mean(delta_i(d_i1==1)==1)]
        frac_b = [frac_b; mean(delta_i(d_i1==1)==0)]

end
end
end

filename = sprintf('output/figures/conditional_alpha=%.1f-gamma=%.3f-delta=%.2f-eta=%.2f.png', alpha, gamma, delta, eta);
saveas(f,filename);
